!Computes the current public spending level gN that makes the government indifferent between repayment and default
! Value A: Repayment (Baseline) with adjusted current gN
! Value B: Default (Baseline)

!INSTRUCTIONS:
!1. To compute policy fcns with (na,nb), need to reset nb2, re-run func first and keep txt files for policy fcns in folder output
!2. Need to call FUNCALLER.F90 beforehand to upload all parameter values

SUBROUTINE WELFARE_GAIN_DEF_YT
!Compute adjusted gN as function of yT

USE GLOBAL
USE OMP_LIB

IMPLICIT NONE

REAL(8),DIMENSION(nb,na)::v_matrixA,v0_matrixA,vaut_matrixA
REAL(8),DIMENSION(nb,na)::b_next_matrixA,q_matrixA,q_matrix_nodefA
REAL(8),DIMENSION(nb,na)::x_matrixA,gn_matrixA,xaut_matrixA,gnaut_matrixA
INTEGER,DIMENSION(nb,na)::default_decisionA

REAL(8),DIMENSION(nb,na)::cn_matrixA,cnaut_matrixA,ct_matrixA,ctaut_matrixA,c_matrixA,caut_matrixA
REAL(8),DIMENSION(nb,na)::tau_matrixA,tauaut_matrixA

REAL(8),DIMENSION(nb,na)::v_matrixB,v0_matrixB,vaut_matrixB
REAL(8)::dev_v,dev_v2,gain_c,criteria
REAL(8)::u_value,a_initial,b_initial,aux0

REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb)
REAL(8)::exp_vA,exp_vB,a00,b00 
INTEGER,PARAMETER::nbA=nb2,nbB=nb2
REAL(8)::distAaccess(nbA,na),distBaccess(nbB,na)
REAL(8)::c0,gn0,slope

INTEGER ILEFT, IRIGHT,NINTV,choice_cond,i_a_global,i_b_global,iter,iter2,nbtemp
EXTERNAL DNORIN
INTEGER::i,i_a,ii_a,i_b,ii_b,iii_b,i_d,i_a0,i_b0,i_a00,i_b00,i_default_today,i_defaultA,convergence,convergence2
REAL(8)::gain_c_new,dev_gain ,valuex

REAL(8)::gain_c_state(na),v0_fun_welfare,weight_gain,v_valor
INTEGER::nbaux,index_opt,i_default_global

REAL(8)::new_value,aahat,aa1,aa0,a01,a02 ,lag_lim 
REAL(8)::valueA,valueB, slope_aux(nb),gn_aus(nb),b_next0
REAL(8)::constant_gn_new,dev_v2_new,constant_gn_old,dev_v2_old 
REAL(8)::objective_fun_gN,q_fun_base,bb
EXTERNAL::objective_fun_gN,q_fun_base
CALL compute_bgrid

open (110, FILE='graphs\b_grid.txt')
DO i=1,nb
    READ(110,*) bgrid(i)
ENDDO

CLOSE(110)
    
 
open (10, FILE='graphs\v.txt')
open (11, FILE='output\gdef.txt')
open (13, FILE='output\gbnext.txt')
open (114, FILE='output\gh.txt')
open (115, FILE='output\ggn.txt')
open (116, FILE='output\ghaut.txt')
open (117, FILE='output\ggnaut.txt')
open (118, FILE='output\gct.txt')
open (119, FILE='output\gctaut.txt')

do i_a = 1,na
do i_b = 1,nb
    READ(114, *) x_matrixA(i_b, i_a)
    READ(115, *) gn_matrixA(i_b, i_a)                     
    READ(118, *) ct_matrixA(i_b, i_a)
    READ(13, *) b_next_matrixA(i_b, i_a) 
ENDDO
ENDDO

do i_b = 1,nb
do i_a = 1,na
    i_a_global = i_a
    i_b_global = i_b

    READ(10, '(F18.10, X, F18.10, X, F18.10)') v_matrixA(i_b, i_a), &
                v0_matrixA(i_b, i_a), vaut_matrixA(i_b, i_a)
    if (vaut_matrixA(i_b, i_a) > v0_matrixA(i_b, i_a)) then
        default_decisionA(i_b, i_a) =2
    else
        default_decisionA(i_b, i_a) =1
    end if
end do
end do

v_matrixB(:,:) = vaut_matrixA(:,:)

do i_a = 1,na
    READ(116, *) valuex
    xaut_matrixA(:, i_a)=valuex
    READ(117, *) valuex
    gnaut_matrixA(:, i_a)=valuex 
    READ(119, *) valuex
    ctaut_matrixA(:, i_a)=valuex 
enddo    

CLOSE(10)
CLOSE(13)
CLOSE(114)
CLOSE(115)
CLOSE(116)
CLOSE(117)
CLOSE(118)
CLOSE(119)

PRINT*,' '
PRINT*,'SOLUTION UPLOADED FOR WELFARE ANALYSIS OF DEF VS. REPAYMENT'
PRINT*,' '

!upload value and policy fcns
PRINT*,'Compute CONDITIONAL consumption gain' 
PRINT*,' '

OPEN(1111,FILE='output//gain_c_state_min_def_yt.txt')

!To apply a_fun, need original value fcns
!Default thresholds computed using model A
vaut_matrix  = vaut_matrixA
v0_matrix    = v0_matrixA
v_matrix     = v_matrixA

index_bfeas4 = 0d0
    
!to compute spline for feasible debt gridpoints
do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt        
ENDDO

index_bfeas2=MAXVAL(index_bfeas4)

do i_a = 1,na         
IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
    ii_b=index_bfeas4(i_a)+1
    nbaux =nb-ii_b+1
    FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	ILEFT  = 0
	IRIGHT = 0
	DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
ENDIF         
ENDDO

do i_a = 1,na
    FDATA = q_matrix_nodef(:,i_a)
	ILEFT  = 0
	IRIGHT = 0
    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
    break_matrix_q(i_a, :) = BREAK_GRID
    coeff_matrix_q(i_a, :,:) = CSCOEF
end do

criteria=1d-5
iter=0
gain_c_state = 0.01d0


a00 = mua/(1d0-rhoa)      !start from unconditional mean
i_a00 = LOCATE2(agrid,a00)

b00 = -1.108d0*(lambda+r)
i_b00 = LOCATE2(bgrid,b00)

!Welfare gains as function of current yT given current debt set to average level
DO ii_a=1,na
    i_a0 = ii_a
    
DO ii_b=i_b00,i_b00
    i_b0 = ii_b

    a_initial = agrid(i_a0)
    b_initial = bgrid(i_b0)
    
    i_default_global=1 !Country does not default

    IF (MIN(v_matrixB(i_b0, i_a0),v0_matrixA(i_b0, i_a0)).LE.-500d0) THEN
        gain_c_state(i_a0) = -1d0
        CYCLE
    ENDIF

    convergence = -1
    convergence2= -1

    gain_c = 0.00d0
    
    !constant_gN is global variable
    constant_gn = gn_matrixA(i_b0, i_a0) 
    dev_gain    = 1d0

    !v_matrixB: Default, baseline
    !v_matrixA: Repayment, adjusted gN
    dev_v2 = v_matrixB(i_b0, i_a0)-v0_matrixA(i_b0, i_a0)
    PRINT*,'Indices of current state (yT,debt) =',i_a0,i_b0
    PRINT*,'Values Repayment, Default=',  v_matrixB(i_b0, i_a0), v0_matrixA(i_b0, i_a0)
    PRINT*,'Difference in Values =', dev_v2
    PRINT*,' '
    
    IF (dev_v2>0) THEN
        gn_aus(ii_a) = 100d0
        PRINT*,'ii_a, govt spending',ii_a, gn_aus(ii_a)
        WRITE(1111,' (F18.10, X, F18.10, X, F18.10, X, F18.10, X, F18.10) ')  agrid(i_a0),  gn_aus(ii_a), gn_matrixA(i_b0, i_a0) ,  v_matrixB(i_b0, i_a0), v_matrixA(i_b0, i_a0)   
    ELSE       

        dev_v2_new = dev_v2
        constant_gn_new = gn_matrixA(i_b0, i_a0) 

        constant_gn = 0.05d0*gn_matrixA(i_b0, i_a0) 
        b_next0= b_next_matrixA(i_b0, i_a0) 
        v_valor=-objective_fun_gN(b_next0,a_initial,b_initial,i_default_global)

        dev_v2_old=v_matrixB(i_b0, i_a0)-v_valor
        constant_gn_old=constant_gn
        
        IF (dev_v2_old<0) THEN
            PRINT*,'problem with initial constant_gn',dev_v2_old
            PRINT*,'vB,vA',v_matrixB(i_b0, i_a0),v_valor
            PRINT*,'i_a0,i_b0',i_a0,i_b0
            PAUSE
            STOP
        ENDIF
        
        iter=0
        
        do WHILE(iter<1000)
                iter=iter+1
                iter2 = 0
                convergence=-1
                iter2 = iter2+1    
                
                !constant_gN is global variable
                constant_gn = (constant_gn_new+constant_gn_old)/2d0

                !use optimize_gN to set current gN while letting transfers be chosen optimally
                call optimize_gN(b_next0, v_valor,a_initial,b_initial,i_default_global)
                dev_v2 =  v_matrixB(i_b0, i_a0) - v_valor
                  
                 IF (dev_v2*dev_v2_old.le.0d0) THEN
                    constant_gn_new=constant_gn
                    dev_v2_new = dev_v2
                ELSE
                    constant_gn_old=constant_gn
                    dev_v2_old = dev_v2
                ENDIF
                
                IF (DABS(constant_gn_new-constant_gn_old)<1d-9) THEN
                    constant_gn=(constant_gn_new+constant_gn_old)/2d0
                    gn_aus(ii_a) = constant_gn
                    PRINT*,'ii_a, govt spending',ii_a, gn_aus(ii_a)
                    WRITE(1111,' (F18.10, X, F18.10, X, F18.10, X, F18.10, X, F18.10) ')  agrid(i_a0),  constant_gn,  gn_matrixA(i_b0, i_a0) ,v_matrixB(i_b0, i_a0), v_matrixA(i_b0, i_a0)   
                    EXIT
                ENDIF

                WRITE(*, '(I5, X, I5, X, F12.8, X, F12.8, X, F12.8, X, F12.8)') ii_b,iter, constant_gn, dev_v2,constant_gn_new,constant_gn_old

    ENDDO   !DO WHILE
    ENDIF
ENDDO       !ii_b
    
ENDDO       !ii_a


CLOSE(1111)

PRINT*,'ADJUSTED GN COMPUTED AS FUNCTION OF YT'


END SUBROUTINE WELFARE_GAIN_DEF_YT
!********************************************************************************************
SUBROUTINE WELFARE_GAIN_DEF_DEBT
!Compute adjusted gN as function of debt

USE GLOBAL
USE OMP_LIB

IMPLICIT NONE

REAL(8),DIMENSION(nb,na)::v_matrixA,v0_matrixA,vaut_matrixA
REAL(8),DIMENSION(nb,na)::b_next_matrixA,q_matrixA,q_matrix_nodefA
REAL(8),DIMENSION(nb,na)::x_matrixA,gn_matrixA,xaut_matrixA,gnaut_matrixA
INTEGER,DIMENSION(nb,na)::default_decisionA

REAL(8),DIMENSION(nb,na)::cn_matrixA,cnaut_matrixA,ct_matrixA,ctaut_matrixA,c_matrixA,caut_matrixA
REAL(8),DIMENSION(nb,na)::tau_matrixA,tauaut_matrixA

REAL(8),DIMENSION(nb,na)::v_matrixB,v0_matrixB,vaut_matrixB
REAL(8)::dev_v,dev_v2,gain_c,criteria
REAL(8)::u_value,a_initial,b_initial,aux0

REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb)
REAL(8)::exp_vA,exp_vB,a00,b00 
INTEGER,PARAMETER::nbA=nb2,nbB=nb2
REAL(8)::distAaccess(nbA,na),distBaccess(nbB,na)
REAL(8)::c0,gn0,slope

INTEGER ILEFT, IRIGHT,NINTV,choice_cond,i_a_global,i_b_global,iter,iter2,nbtemp
EXTERNAL DNORIN
INTEGER::i,i_a,ii_a,i_b,ii_b,iii_b,i_d,i_a0,i_b0,i_a00,i_b00,i_default_today,i_defaultA,convergence,convergence2
REAL(8)::gain_c_new,dev_gain ,valuex

REAL(8)::gain_c_state(na),v0_fun_welfare,weight_gain,v_valor
INTEGER::nbaux,index_opt,i_default_global

REAL(8)::new_value,aahat,aa1,aa0,a01,a02 ,lag_lim 
REAL(8)::valueA,valueB, slope_aux(nb),gn_aus(nb),b_next0
REAL(8)::constant_gn_new,dev_v2_new,constant_gn_old,dev_v2_old 
REAL(8)::objective_fun_gN,q_fun_base,bb
EXTERNAL::objective_fun_gN,q_fun_base
CALL compute_bgrid

open (110, FILE='graphs\b_grid.txt')
DO i=1,nb
    READ(110,*) bgrid(i)
ENDDO

CLOSE(110)
 
open (10, FILE='graphs\v.txt')
open (11, FILE='output\gdef.txt')
open (13, FILE='output\gbnext.txt')
open (114, FILE='output\gh.txt')
open (115, FILE='output\ggn.txt')
open (116, FILE='output\ghaut.txt')
open (117, FILE='output\ggnaut.txt')
open (118, FILE='output\gct.txt')
open (119, FILE='output\gctaut.txt')

do i_a = 1,na
do i_b = 1,nb
    READ(114, *) x_matrixA(i_b, i_a)
    READ(115, *) gn_matrixA(i_b, i_a)                     
    READ(118, *) ct_matrixA(i_b, i_a)
    READ(13, *) b_next_matrixA(i_b, i_a) 
ENDDO
ENDDO

do i_b = 1,nb
do i_a = 1,na
    i_a_global = i_a
    i_b_global = i_b

    READ(10, '(F18.10, X, F18.10, X, F18.10)') v_matrixA(i_b, i_a), &
                v0_matrixA(i_b, i_a), vaut_matrixA(i_b, i_a)
    if (vaut_matrixA(i_b, i_a) > v0_matrixA(i_b, i_a)) then
        default_decisionA(i_b, i_a) =2
    else
        default_decisionA(i_b, i_a) =1
    end if
end do
end do

v_matrixB(:,:) = vaut_matrixA(:,:)

do i_a = 1,na
    READ(116, *) valuex
    xaut_matrixA(:, i_a)=valuex
    READ(117, *) valuex
    gnaut_matrixA(:, i_a)=valuex 
    READ(119, *) valuex
    ctaut_matrixA(:, i_a)=valuex 
enddo    

CLOSE(10)
CLOSE(13)
CLOSE(114)
CLOSE(115)
CLOSE(116)
CLOSE(117)
CLOSE(118)
CLOSE(119)

PRINT*,' '
PRINT*,'SOLUTION UPLOADED FOR WELFARE ANALYSIS OF DEF VS. REPAYMENT'
PRINT*,' '

!upload value and policy fcns
PRINT*,'Compute CONDITIONAL consumption gain' 
PRINT*,' '

OPEN(1111,FILE='output//gain_c_state_min_def_debt.txt')

!To apply a_fun, need original value fcns
!Default thresholds computed using model A
vaut_matrix  = vaut_matrixA
v0_matrix    = v0_matrixA
v_matrix     = v_matrixA

index_bfeas4 = 0d0
    
!to compute spline for feasible debt gridpoints
do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt        
ENDDO

index_bfeas2=MAXVAL(index_bfeas4)

do i_a = 1,na         
IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
    ii_b=index_bfeas4(i_a)+1
    nbaux =nb-ii_b+1
    FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	ILEFT  = 0
	IRIGHT = 0
	DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
ENDIF         
ENDDO

do i_a = 1,na
    FDATA = q_matrix_nodef(:,i_a)
	ILEFT  = 0
	IRIGHT = 0
    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
    break_matrix_q(i_a, :) = BREAK_GRID
    coeff_matrix_q(i_a, :,:) = CSCOEF
end do

criteria=1d-5
iter=0
gain_c_state = 0.01d0


a00 = mua/(1d0-rhoa)      !start from unconditional mean
i_a00 = LOCATE2(agrid,a00)

b00 = -1.108d0*(lambda+r)
i_b00 = LOCATE2(bgrid,b00)
    
!Welfare gains as function of current debt given yT set to unconditional mean
DO ii_a=i_a00,i_a00   
    i_a0 = ii_a
    
DO ii_b=1,nb
    i_b0 = ii_b

    a_initial = agrid(i_a0)
    b_initial = bgrid(i_b0)
    
    i_default_global=1 !Country does not default

    IF (MIN(v_matrixB(i_b0, i_a0),v0_matrixA(i_b0, i_a0)).LE.-500d0) THEN
        gain_c_state(i_a0) = -1d0
        CYCLE
    ENDIF

    convergence = -1
    convergence2= -1

    gain_c = 0.00d0
    
    !constant_gN is global variable
    constant_gn = gn_matrixA(i_b0, i_a0) 
    dev_gain    = 1d0

    !v_matrixB: Default, baseline
    !v_matrixA: Repayment, adjusted gN
    dev_v2 = v_matrixB(i_b0, i_a0)-v0_matrixA(i_b0, i_a0)
    PRINT*,'Indices of current state (yT,debt) =',i_a0,i_b0
    PRINT*,'Values Repayment, Default=',  v_matrixB(i_b0, i_a0), v0_matrixA(i_b0, i_a0)
    PRINT*,'Difference in Values =', dev_v2
    PRINT*,' '
    
    IF (dev_v2>0) THEN
        gn_aus(ii_a) = 100d0
        PRINT*,'ii_a, govt spending',ii_a, gn_aus(ii_a)
        WRITE(1111,' (F18.10, X, F18.10, X, F18.10, X, F18.10, X, F18.10) ')  bgrid(i_b0),  gn_aus(ii_a), gn_matrixA(i_b0, i_a0) ,  v_matrixB(i_b0, i_a0), v_matrixA(i_b0, i_a0)   
    ELSE       

        dev_v2_new = dev_v2
        constant_gn_new = gn_matrixA(i_b0, i_a0) 

        constant_gn = 0.05d0*gn_matrixA(i_b0, i_a0) 
        b_next0= b_next_matrixA(i_b0, i_a0) 
        v_valor=-objective_fun_gN(b_next0,a_initial,b_initial,i_default_global)

        dev_v2_old=v_matrixB(i_b0, i_a0)-v_valor
        constant_gn_old=constant_gn
        
        IF (dev_v2_old<0) THEN
            PRINT*,'problem with initial constant_gn',dev_v2_old
            PRINT*,'vB,vA',v_matrixB(i_b0, i_a0),v_valor
            PRINT*,'i_a0,i_b0',i_a0,i_b0
            PAUSE
            STOP
        ENDIF
        
        iter=0
        
        do WHILE(iter<1000)
                iter=iter+1
                iter2 = 0
                convergence=-1
                iter2 = iter2+1    
                
                !constant_gN is global variable
                constant_gn = (constant_gn_new+constant_gn_old)/2d0

                !use optimize_gN to set current gN while letting transfers be chosen optimally
                call optimize_gN(b_next0, v_valor,a_initial,b_initial,i_default_global)
                dev_v2 =  v_matrixB(i_b0, i_a0) - v_valor
                  
                 IF (dev_v2*dev_v2_old.le.0d0) THEN
                    constant_gn_new=constant_gn
                    dev_v2_new = dev_v2
                ELSE
                    constant_gn_old=constant_gn
                    dev_v2_old = dev_v2
                ENDIF
                
                IF (DABS(constant_gn_new-constant_gn_old)<1d-9) THEN
                    constant_gn=(constant_gn_new+constant_gn_old)/2d0
                    gn_aus(ii_a) = constant_gn
                    PRINT*,'ii_a, govt spending',ii_a, gn_aus(ii_a)
                    WRITE(1111,' (F18.10, X, F18.10, X, F18.10, X, F18.10, X, F18.10) ')  bgrid(i_b0), constant_gn,  gn_matrixA(i_b0, i_a0) ,v_matrixB(i_b0, i_a0), v_matrixA(i_b0, i_a0)   
                    EXIT
                ENDIF

                WRITE(*, '(I5, X, I5, X, F12.8, X, F12.8, X, F12.8, X, F12.8)') ii_b,iter, constant_gn, dev_v2,constant_gn_new,constant_gn_old

    ENDDO   !DO WHILE
    ENDIF
ENDDO       !ii_b
    
ENDDO       !ii_a


CLOSE(1111)


PRINT*,'ADJUSTED GN COMPUTED AS FUNCTION OF DEBT'

END SUBROUTINE WELFARE_GAIN_DEF_DEBT



